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Abstract 

We describe a general analytic-numerical reduction scheme for evaluating 
any two-loop diagrams with general kinematics and general renormalizable in- 
teractions, whereby ten special functions form a complete set after tensor reduc- 
tion. We discuss the symmetrical analytic structure of these special functions 
in their integral representation, which allows for optimized numerical integra- 
tion. The process Z — » bb is used for illustration, for which we evaluate all the 
three-point, non-factorizable g 2 a s mixed electroweak-QCD graphs, which de- 
pend on the top quark mass. The isolation of infrared singularities is detailed, 
and numerical results are given for all two-loop three-point graphs involved in 
this process. 
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1 Introduction 

The success of the perturbative aspect of the Standard Model is truly impressive, 
the theory being now tested at higher than one-loop accuracy. It is natural that 
over the past years various methods have been proposed to calculate higher loop 
graphs. Beyond one loop, there are situations in which one can neglect internal 
masses, or exploit the kinematics to evaluate amplitudes analytically at some special 
point where the Feynman diagrams become simpler. These special situations, while 
ingenious, cover only certain classes of physical processes. 

There have been attempts to formulate solutions for two-loop Feynman diagrams, 
which ideally should be as general as the one-loop solutions are. For the massless case 
a lot of progress has been reached recently Jl|. On the other hand, for the general 
massive case, it has become clear that some numerical treatment appears unavoidable 
because of the complexity of the scalar integrals involved. 
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In refs. [0, |3| a reduction scheme was found for treating massive self-energy dia- 
grams at two-loop, with the resulting master scalar integrals being evaluated numer- 
ically. 

In ref. JiJ we investigated the problem of a general massive two-loop algorithm, 
which would deal with multi-leg diagrams as well. Once we subscribe to a semi- 
numerical approach, some demands must be made, which are based on the generality, 
the universality, the effectiveness and the accuracy that a particular formalism engen- 
ders. We have shown that for any external kinematics and internal masses, we can 
reduce every two-loop amplitude for a process due to any renormalizable interactions 
into a set of ten scalar integrals and their derivatives. 

It is not hard to see how a result like this is obtained if we consider first a scalar 
theory with trivial interactions |J. Then, there would be two internal momenta, p 
and q, and by the use of Feynman parameters to combine sets of denominators in 
which either p or q or p + q runs through, and to make shifts in p and q, we have an 
integrand proportional to 



[(p + k) 2 + mj ] Ql (q 2 + m|)° 2 [(p + q) 2 + m|] a3 

in which a four-vector k is left over, which is a linear function of the external momenta 
and Feynman parameters. The "masses" m\ 2 3 are functions of external momenta, 
internal masses and Feynman parameters. These are to be integrated over p, q and 
a set of Feynman parameters. By partial differentiating with respect to m 2 2 3 , there 
is in principle only one basic function arising from a\ = a 2 = 012 = 1 we need to 
know, although for convenience one may add a\ — 2, a 2 — 0:3 = 1. It is clear that 
different interactions and different graphs will give different polynomials of Feynman 
parameters to the numerators, and also different k and m 2 2 3 . An extension of this 
construction to tensor integrals will give us the main result mentioned earlier. 

In the following section, we discuss the steps and arguments needed to perform a 
tensorial decomposition. We show that there is a small set of basic functions which are 
needed, which have simple, one- dimensional integral representations. Their analytic 
structure is easy to see, and therefore the integration contour can be extended into 
the complex plane and optimized for rapid numerical convergence. In Section 3, we 
shall give a detailed exposition of the numerical techniques involved. 

Given that two-point functions are usually simpler to calculate, many calcula- 
tions existing in the literature used unitarity cuts of self-energies to calculate lower 
loop-order inclusive decay rates. We would like to stress that, given a process, our 
formalism allows us to calculate the amplitudes individually, as opposed to using uni- 
tarity cuts. As an immediate consequence, it becomes possible to obtain rates for 
exclusive processes, without the need to integrate over the whole final-state phase 
space. 

In section 4, for illustration, we discuss in detail the application of our method to 
the important physical process Z — > bb. 
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Figure 1: Expressing generic massive two- loop Feynman diagrams as integrals over 
sunset-type functions. 

2 Analytical reduction 

2.1 Relation with sunset-type integrals 

The starting point of calculating a massive two-loop diagram is to express it in terms 
of integrals of a standard type. Topologically, any generic two-loop diagram can be 
represented in the form shown in figure 1. By introducing a set of Feynman parameters 
X to combine all propagators which have the same loop integration momentum p, q, 
or r = p + q, the graph can always be represented as an integral over a sunset-type 
two-loop integral. This is illustrated in figure [3]. All dependence on the external 
momenta fcj and internal masses rrii is now contained in the variables mf, m|, mf, 
and k, which also depend on the set of Feynman parameters X. 

The original Feynman diagram is written at this point as an integral over tensor 
functions of the following type (all momenta are rotated into Euclidean): 



By casting the graph in this form, our strategy is to develop a uniform treatment of 
all possible sunset-type functions which can be generated from generic two-loop Feyn- 
man graphs. The further tensor reduction and decomposition into standard scalar 
functions is done by using standard formulae common for all diagrams. Therefore it 
can be automatized in the form of an algebraic manipulation program. 

At the end, the remaining integral over the Feynman parameters X, represented 
in figure |l], is performed numerically. 




(2) 
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2.2 Tensor reduction 



Tensor integrals of the type in eqn. § need to be decomposed into scalar integrals. By 
Lorentz covariance, this two- loop integral is a tensor constructed from the external 
momentum and the metric tensor g^ v . Given that, one way of obtaining the tensor 
decomposition would be to write down all tensor structures allowed by the symmetry 
of the integral, use appropriate projectors, and solve the resulting equations. 

Another way to do this is by decomposing the loop momenta p and q into com- 
ponents parallel and orthogonal to the external momentum 

p£=p"-^*" , 4L = <f-^w • (3) 

After this decomposition, the tensor decomposition is obtained by noticing that the 
functions with an odd number of transverse loop momenta p^_ and vanish, while 
the even functions are transverse to k^ [|J. In ref. j|] we have shown that the resulting 
scalar coefficients of the tensor decomposition are integrals of the following form: 



K b ia2a3 (m u m 2 ,m 3 ;k 2 ) = Jd n pd n q- 



(p-k) a (q-kY 



[(p + k) 2 + m 2 ]^ (q 2 + m 2 .) 02 (r 2 + m|) a s 

(4) 

We give in the following the tensor decompositions for all tensor integrals of the 
type of eq. || of rank 1,2, and 3: 



[211] 1 1 ' [211] z 1 ' [211] 

= t»\A x + gTiA* , p^j = tV b Ai + g»\A % , |^ = t»\A x + g»\A 2 



q^p v p x 
p^q v q x 



q^q v q x 



(T^k x + r^k" + T uX k") 7 A 1 + {g^k x + g» x k v + g vX k%A 2 

(r^k x + r^k" + T uX k") s A 1 + (g» v k x + g» x k v + g uX k») 8 A 2 

(T^k x + r^k" + T uX k") 9 A 1 + (g^ u k x + g» x k v + g vX k») 9 A 2 

+ ( g ^k x + g^ ~2g^) 9 A, 
(r^k x + r» x k v + T^k^ioA! + (g^k x + g» x k u + g uX k^) w A 2 



Tensor decompositions for more than three Lorentz indices are derived in a similar 
way. In the formulae above, a loop integration J d n p d n q is understood, and we used 
the following notations: 
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[211] = i(p + k) + m\} 2 (q 2 + m\) [r 2 + m 2 ) 



(5) 



and: 



1A1 
aA\ 

5^1 
6^1 
7^1 
8^1 
9^1 
9^3 
10^1 



Hi 



1 




n 




k 2 


n 




1 


1 




n 




A; 2 


n 




1 


1 




n 




k 2 


n 




1 



7^4 

n + 2 
fcV 3(n- 1) 
n + 2 



A: 2 / 3(n-l) 
1 \ 2 n + 2 



1 1 



A; 2 / 3(n-l) 

' 211 T / 211 



7*7 



n 



n 



2^1 

4 A 2 : 

5^2 ; 
6^2 : 
7^2 -- 

& A 2 
9 A 2 - 
■ (H 5 + Hq 



¥ H2 

k 2 I 211 
1 

k 2 
1 



k 2 ' 211 



2 



k 2 ' 211 
1 \ 2 1^ 



k 2 ) 3' 211 

1-621 

A; 2 / 3 /m 

1 \ 2 1 - 
— I --P 12 
A; 2 ' ' 



211 



1 



n + 2 



A; 2 / 3(ra-l) 



10 



1(A= ij_yip«3 



211 



A; 2 



7i 3 



In the expressions above, we introduced a set of ten scalar integrals Hi, which are 
related to the V$u functions, and are defined as follows: 



n 2 
n 3 

7*4 

7*5 

7*6 



d n pd n q 
J d n pd n q 
J d n pd n q 
d n pd n q 
J d n pd n q 
d n pd n q 



1 



[(p + A;) 2 + m 2 } 2 (q 2 + m|) (r 2 + m|) 
p ■ k 

[(p + A;) 2 + m 2 ] 2 (q 2 + m 2 ,) (r 2 + m|) 
q ■ k 

[(p + A;) 2 + m 2 ] 2 (q 2 + m 2 ) (r 2 + ml) 

(p ■ k) 2 — ^Aj 2 p 2 
[(p + A;) 2 + m 2 ] 2 (q 2 + m 2 ,) (r 2 + ml) 

(p-k)(q-k)-±k 2 (q-p) 
[(p + A;) 2 + m 2 ] 2 (g 2 + m|) (r 2 + mi) 

(g ■ A;) 2 - \k 2 q 2 

[5 + A;) 2 + m 2 ] 2 (g 2 + ml) (r 2 + mi) 



5 



• k f - ^k 2 p 2 (p ■ k) 



l~L-j = / d n p d n q 

J [(p + k) 2 + m 2 ] 2 (q 2 + m 2 ,) (r 2 + m 2 ) 

Hd = J dydnq {p-k){q-k) 2 -^- 2 k 2 [2{p.q){ q .k) + q 2 {p.k)} 



(p-k) 2 (q-k)-^k 2 p 2 (q-k) 
[ip + k) 2 + m 2 ] 2 (q 2 + m 2 ) (r 2 + m|) 



[{p + k) 2 + m 2 ] 2 (g 2 + m 2 .) (r 2 + m 2 ) 

f n n (q- k) 3 - -^k 2 q 2 {q ■ k) 

Hio = J d n pd n q jj- - fc)2 - m 2] 2 ( ?2 + m 2) p + m 2) ( 6 ) 

As it can be seen in eqs. 5, these scalar functions appear naturally in the tenso- 
rial decomposition. As discussed in the following section, the functions Tii are only 
logarithmically divergent in the ultraviolet. Because of this, they have very simple 
integral representations which can be used for numerical computation. Once a way 
of calculating the special functions Tii is available, the V% \ x functions can easily be 
recovered by partial fractioning eqns. 6. The partial fractioning generates essentially 
trivial products of one-loop tadpoles. The conversion formulae are given in Appendix 
A. 



2.3 Integral representation of the Hi functions 

Because the ultraviolet behaviour of the functions Tii is logarithmic, fairly simple and 
symmetric integral representations can be found [ffl: 



n 8 



71 



2 1. 1 7T 2 2 

^ - -I 1 - 2 7mi) - 2 + ^2 " 7mi + 7m i + 1 



7T A k 2 
7rV 



2 11 

"^ + 7 ( 2 _27mi) + 



13 _ 7T 2 

~8 ~ 12 



+ 



7m! 



7mj - h 1 



1 1/1 \ 13 7T 2 J mi 



e v 4 



4/j.2\2 



7T 4 (fc 2 ) 2 



Tl\k 2 ) 2 
TT\k 2 f 
TT\k 2 f 



16 24 

_2 



4 



3 13 7mi _175 vr 2 3 7 2 3 
2e 2 + 6 2 96 + 16 + 4 + 4^ 



A _ 1 3 ^i 175 7T 2 3 7 ^ 3 
4e 2 e 4 192 32 8 4 5 



1 1/1 7mi 

2? ~ e 24 - ~2~ 

1 15 
~^"7 ( 24 +7mi 



2e 2 + eM8 + 



19 


7T 2 


7m 1 | 

24 


32' 


"48" 


287 


7T 2 


57mi 


192 


24 


24 


287 
384 


7T 2 

+ 48 


57 mi 
48 



y 2 3 

4 4 6 



+ 



2 7 



4 2 fc 



6 



ir\k 2 f 
n 4 (k 2 ) 3 



1 

3? 



1 / 1 7mi x 95 



— -(— 

4e 2 + e^96 



e v 24 3 



192 



7T 

72 



283 
768 



Tm i 

24 

Tmi 



7 2 



1 



7T 

96 4 9G 



7 2 

/ mi 



6 2 
1 
2 



(7) 



Here, n = 4 + e is the space-time dimension, and 7 m = 7 + log(7rm 2 //z 2 ) (7 is 
the Euler constant, and /ii is the 't Hooft mass). The special functions hi which 
appear in the formulae above are the finite part in the 1/e expansion of Tii. They are 
defined by the following integral representations. Except for special values of their 
arguments, they cannot be further integrated into well-studied functions, such as the 
familiar polylogarithms, and as such, our strategy is to evaluate them numerically 
directly from their integral representations: 



hi(m!, 


m 2 , 


m 3 ; 


k 2 ) = 


/ dxg(x) 
Jo 






h 2 (m 1 , 


m 2 , 


m 3 ; 


k 2 ) = 


/ dx \g(x) 
Jo 


+ fi(x)] 




h 3 (m 1 , 


m 2 , 


m 3 ; 


k 2 ) = 


/ dx [g(x) 
Jo 


+ f 1 (x)](l-x) 




^K, 


m 2 , 


m 3 ; 


k 2 ) = 


/ dx \g(x) 
Jo 


+ fi(x) + f 2 (x)\ 




/i 5 (mi, 


m 2 , 


m 3 ; 


k 2 ) = 


/ dx \g(x) 
Jo 


+ h(x) + f 2 (x)](l-x) 




h(mi, 


m 2 . 


m 3 ; 


k 2 ) = 


/ dx [g(x) 
Jo 


+ h(x) + f 2 (x)](l-x) 2 




h 7 (nii, 


m 2 . 


m 3 ; 


k 2 ) = 


/ dx \g(x) 
Jo 


+ fi(x) + f 2 (x) + f 3 (x)} 




hsirri!, 


m 2 . 


m 3 ; 


k 2 ) = 


/ dx [g(x) 
Jo 


+ f 1 (x) + f 2 (x) + f 3 (x)](l 


-x) 


hg(™>i, 


m 2 , 


m 3 ; 


k 2 ) = 


/ dx \g(x) 
Jo 


+ Mx) + f 2 (x) + f 3 (x)](l 


-xf 


/iio(mi, 


m 2 . 


m 3 ; 


k 2 ) = 


I dx \g(x) 
Jo 


+ A(x) + f 2 (x) + f 3 (x)](l 


-xf 



The four building blocks g(x), fi(x), f 2 (x), and f 3 (x) of these one-dimensional 
integral representations are given in Appendix B. 

2.4 Completeness of the {Hij^ijQ special functions for renor- 
malizable theories 

The ten special functions 7ii we introduced in the previous section, or, equivalently, 
their corresponding V^ b lCt2a , are sufficient for treating all two- loop Feynman graphs 
which may be encountered in renormalizable theories. 
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To see this, it is useful to define the following auxiliary two-loop functions: 



■V ab <m m m ■ h 2 \ - fj n -nJ n n (j? • K) (g • fcj 

I aia2a3 {m u m 2 ,m 3 ,K )- J a pa ? (p2 + m 2 )ai (g2 + m 2 )Q2 [(r + fc)2 + m 2 ]Q3 " 

(9) 

Vaia 2 a 3 are trivially related to V^ b ia2a3 by a simple redefinition of the loop momen- 
tum (see Appendix A). The only reason for introducing them is for simplifying the 
discussion of this section; some recursion relations can be written more compactly in 
terms of these scalar functions. 

To prove our assertion that the set of ten {Hi} functions is sufficient for the case of 
renormalizable theories, we notice that not all functions V^ ia2a of various indices are 
independent. There are recursion relations which relate functions of different indices. 

We will define the "degree" of a V^_ a a function as a\ + a 2 + «3 — a — b. The 
degree can be increased by differentiating with respect to the mass variables: 

1 d 

K^+ia 2 a 3 (mum 2 ,m 3 ;k 2 ) = - — ^i 2a ,(mi,m 2) m 3 ; k 2 ) , (10) 

and similarly for a 2 and a 3 . 

Functions of the same degree are related by recursion relations obtained by dif- 
ferentiating with respect to the external momentum variable k^\ 



-pa+l b 

ai+1 «2 «3 



-T^a 6+1 



1 



2 ai 
1 

2^~ 2 



2k z 



2k 



d 
dk 2 
,_d_ 

dk 2 



(a + b) 
(a + b) 



nk 2 

-pab , Z2_-pa-lb 

Ct\ «2 Q3 2q/ Ql 012 ° 3 

hk 2 

<pab 1 Ufv <pab-l 

Ot\ OL2 O-Z ' a l Q 2 03 ' 



and 



2k 



_d_ 

dk 2 



(a + b) 



V, 



a b 

ai a.2 03 



-2a 3 



k 2 V, 



a b 

Ql Q2 03- 



1 pa+l 6 1 <p< 

-1 1 ' ai ai «3+l 1 aict2«3+l 



■>ao+l 



(12) 

When calculating a specific process, it is thus possible to just focus on the low- 
est order V^ b a20l3 functions involved. All the other can be derived from these by 
differentiation, by using the relations above. 

At the same time, in renormalizable theories, there is a lower bound on the possible 
degree of V^ a2a3 functions involved in two-loop graphs, imposed by the dimension 
of the operators in the Lagrangian. The minimal degree is one, and is attained for 
instance by a two-loop diagram such as the one shown in figure ^. Additional external 
legs can at most increase the degree of the V^ b a2a3 functions involved. 

Therefore, the most obvious choice of a basic set is the V^ b a2a:j functions with 
ol\ — 0L2 — &3 — 1) ^ + b — 0, 1, 2: 



8 



Figure 2: Two-loop diagram which involves V^ a2a functions of minimal degree one 
It consists of a fermion loop, a boson propagator, and non- derivative vertices. 



degree = 3 
degree = 2 
degree = 1 



poo 
' in 
•pio 
' in 

"P2 
' 111 



poi 

' ill 
pll 
' 111 



(13) 



1 111 



However, we prefer to use the functions with ol\ 
instead: 



2, a 2 = a 3 = 1, a + b = 0, 1, 2, 3 



degree = 4 
degree = 3 
degree = 2 
degree = 1 



poo 

' 211 
plO 
' 211 
p2 
' 211 
p3 
' 211 



pOl 



211 



pll 



211 



p0 2 



(14) 



211 



p21 



211 



pi 2 



211 



p0 3 



211 



It turns out that this equivalent set of functions has simpler integral representa- 
tions than the functions of eqs. |13| When needed, the eqs. [13] functions can be derived 
from eqs. [T4] by partial p: 



V, 



a b 

ai OL2 «3 



1 



n - (ai\ + a 2 + a 3 ) + (a + b)/2 



Ol2 a-j, + l 



2k 2 



d_ 

dk 2 



{a + b) 



<pab 

a 1 ct2 «3 



(15) 



Because the set of functions of eqs. [b| and the ten Hi functions differ essentially by 
simple, one-loop tadpole contributions given in Appendix A, this proves our assertion 
that the ten 7ij functions are sufficient for treating renormalizable theories at two- 
loop. 

We would like to stress that additional recursion relations and symmetries exist 
among this set of functions. An example of such a symmetry relation is given in eq. 
16 of ref. @], and additional ones may be present. Such relations can be used, in 
principle, for restraining the number of master functions involved. However, the set 
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of ten Hi that we proposed has the advantage of having clean, symmetrical integral 
representations, where permutations of mass arguments are not involved, while sym- 
metry relations obtained by loop momenta redefinitions often interchange the mass 
arguments. This simplicity is an advantage in practical calculations done by com- 
puter algebra, where a simpler algorithm may be preferable to a more complicated 
one which produces somewhat more compact results. 

Thus, on general grounds, all tensor integrals V® b ia2a3 with a + b > 4 can be 
expressed by the standard set of ten functions via mass or momentum differentiation. 
The differentiations required can be carried out either by hand, or automatically, by 
using a computer algebra program such as Mathematica or Maple. 

However, in practical calculations of Feynman graphs, before performing the ten- 
sor decomposition, it is advantageous to exploit all possible partial fractioning of the 
loop momenta, instead of resorting to recursion relations after the tensor decomposi- 
tion. This can result in more compact results. This is in particular the case with the 
examples given in section 4. There, the calculation can be simplified by noticing that 
in the Feynman gauge, loop momenta p and q in the numerators come from fermion 
propagators, and there are no more than four of them. For those terms with four 
powers of p and q, one can algebraically rearrange them, so that two of these powers 
will contract as p 2 , q 2 or p ■ q, and then perform partial fractioning of these terms, to 
reduce them into P functions with a + b < 3. 

3 Numerical integration 

Once the tensor reduction is done analytically, the original Feynman diagram is ex- 
pressed, according to figure [I], as an integral over the set of remaining Feynman 
parameters X. The integrand itself consists of a sum of hi special functions and 
possibly trivial functions such as logarithms and rational functions of the kinematic 
variables ml, m 2 ,, m 2 , and k 2 . 

The hi functions themselves are given by the one-dimensional integral representa- 
tions given in eqns. |[ Within our method, all these integration steps are performed 
numerically in general. 

Therefore, it is natural to separate the numerical integration into two distinct 
steps: the routines to calculate the hi functions from their integral representations, 
and the final integration over the remaining Feynman parameters X. 

When performing the numerical integration, special care is needed to make sure 
that the integration is performed on the physical sheet. Let us start with the inte- 
gration of the hi functions. 

In figure ^ we plot the functions g(x) and fi(x) as a function of the integration 
parameter x. This is a typical behaviour above the threshold. At x = there is 
an integrable singularity of the logarithmic type, and similarly at x — 1. By mass 
or momentum differentiation the integrable nature of these singularities in x is not 
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/a (a) V 



y Im 



Figure 3: Typical behaviour of functions g(x), fi(x), fe(x), fs(x) and of their deriva- 
tives as a function of the integration parameter x. The specific kinematic variables 
in these plots are above threshold, —k 2 > (mi + m<i + m3) 2 , where an imaginary part 
is present. 

changed; this makes possible the treatment of more complicated topologies with ad- 
ditional external legs by using the same numerical approach. Along the integration 
path there are two branching points, between which the function acquires an imagi- 
nary part. The integration path must be chosen such as to reproduce correctly this 
imaginary part. 

By continuing the integration parameter x into the complex plane, we obtain 
a picture of the hi function's integrand as shown in figure |j. It is clear that if 
the position of the two branching points is known, a correct integration path can be 
calculated automatically by the computer program. Finding the singularities involves 
some subtleties. 

Adopting the notations introduced in Appendix B, and inspecting the expressions 
of the functions g(x), fi(x), fi{x) and fz{x) given there, the branching points we are 
looking for must be among the solutions of the equation: 

A = (l + /t 2 - / u 2 ) 2 + 4K 2 /i 2 -4m 2 r? = . (16) 
There are four solutions: 

£1,2 = [-a + b + nl ± ^ (a - b - /if) 2 - 46/if ] 

x 3,4 = —^[-a + b + /j,l±yj(a-b- /i 2 .) 2 - 4fyt| ] , 

£ 2 = l-/t 2 T 2 v /Z ^ 2 • (17) 
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Real 



Imaginary 




Figure 4: Typical behaviour of the functions g(x), fi(x), foix), f?,(x) and of their 
derivatives in the complex integration parameter x. The integration proceeds from 
x — to x — 1. 

To establish which two of the four solutions are the branching points we are looking 
for, it is useful to note that the causality of the Green's functions can be expressed 
in at least two equivalent ways, which ought to lead to the same prescription for 
the integration path. The causality condition is expressed by the ir] prescription in 
the Feynman propagator, which means to shift all masses in the propagators m 2 — > 
m 2 — irj. An equivalent way to impose causality is to calculate the Euclidian Green's 
functions and to go afterwards to physical momenta, approaching the cut on the 
positive real axis from above. This amounts to shifting the external momentum 
k 2 — > k 2 + if]. These two prescriptions ought to be equivalent, and therefore have to 
fix the location of the physical singularities with respect to the real axis in the same 
way. For x\ and x 2 , at — k 2 > 1 both prescriptions lead to the same change, and 
therefore these are the singularities of the g and fa functions. For and X4, the two 
prescriptions lead to opposite changes in the imaginary direction. Since causality fixes 
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Figure 5: A typical shape for the complex integration path, defined in terms of higher 
order spline functions. Both the path and the Jacobian are smooth functions. 

the location of the singularities of the Green's function uniquely, £3 and £4 cannot 
correspond to real singularities of g and fa. Therefore g and are analytical at these 
two points. The g and /j functions themselves have only two branching points at 
X\ and x 2 , because the singularities at x 3 and x 4 are compensating among the four 
terms of g in eq. 26. It is interesting to notice that the spurious solutions x 3 and x 4 
correspond to the spurious, unphysical thresholds discussed in ref. in the context 
of the relation of the scalar sunset diagram with a generalization of hypergeometric 
functions. 

Once the positions of these two singularities are known, the computer program 
can automatically compute an integration path which avoids them. Because we are 
interested in a high accuracy and efficiency routine, we used an adaptive deterministic 
integration algorithm. Such integration routines are very accurate provided that the 
integrand is smooth enough. The integrand itself is, of course, an analytic function 
along the complex integration path, and to preserve its smoothness it is advantageous 
to define a smooth integration path as well. We use an integration path defined in 
terms of higher-order spline functions such that both the path and its Jacobian are 
smooth functions. A typical path is shown in figure [|. 

Along these lines, a fully automatic computer program can be written, which first 
identifies the singularities, then calculates a suitable complex integration path, and 
then performs the numerical integration to obtain the hi functions starting from their 
integral representations given in eqs. §. By using adaptive integration routines, one 
typical evaluation of an hi function or of a derivative at eight digits takes of the order 
of 30 ms on a 600 MHz Pentium processor. 

In practical calculations of Feynman graphs, a number of mass or momentum 
derivatives of the functions hi become necessary. They can easily be obtained either 
by hand, or automatically, by computer algebra programs such as Mathematica or 
Maple. 

The numerical integration over the remaining Feynman parameters X is carried 
out along similar lines || |8|, |9|. The dimensionality of this final integration depends 
on the topology of the diagram and on the number of legs. 
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Figure 6: Plots of the ten special functions ^(mj, m|, m|; —k 2 ) as a function of the 
external momentum variable k 2 . The plots given here are for m 2 = m\ = m 2 = 1. 

With increasing complexity of the diagram, the methods we describe in this paper 
will still be applicable in principle, but in practice will result in larger integrands 
and a final numerical integration of higher dimension. Thus its applicability will be 
limited in practice by the available computing power and by the ability of handling 
potentially large expressions which result from the tensor reduction in an error free 
way. 



4 Examples 

4.1 Three-point two-loop diagrams contributing to Z — > bb 

Here we would like to illustrate by means of a concrete example how the algorithm 
described above works in practice. Examples involving two-point functions were given 
in ref. 0, [J. Here, we treat all two- loop three-point diagrams which contribute to 
the important physical process Z — > bb at 0(a s g 2 ). The diagrams involved are shown 
in figure [5]. In the complete calculation of this process, there are also self-energy 
type diagrams which contribute to the b wave function renormalization. These two- 
point function graphs are simpler than the three-point graphs both analytically and 
numerically. Within our method, they were calculated in ref. || . The process Z — > bb 
at 0(a s g 2 ) at 0{a s g 2 ) was calculated by means of a mass expansion in ref. J7J. 
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Figure 7: Two-loop three-point diagrams contributing to Z — > bb at 0(a s g 2 ). 
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4.2 Isolating the IR singularities 



The two-loop methods we describe in this paper are intended primarily for the cal- 
culation of massive diagrams. However, sometimes massless particles are involved in 
a calculation and may lead to infrared singularities. We would like to stress that for 
purely massless calculations, such as QCD radiative corrections, more efficient meth- 
ods are available already in the literature, which make use of the absence of masses 
|L|. It is often possible in the massless case to carry out the calculation completely 
analytically. The discussion in this section applies to cases where the main difficulty 
is related to the presence of several masses, while the infrared structure is relatively 
simple. Typical examples are the purely electroweak and the mixed electroweak-QCD 
radiative corrections, such as the Z — > bb at 0(a s g 2 ) process discussed here. 

The diagrams V2, V4, V$, and V? are infrared divergent. The most common treat- 
ment of IR divergencies uses dimensional regularization to separate both the UV and 
the IR singularities. In our approach however, this is not directly possible because in 
the expressions of eqs. 7 the 1/e expansion already separated the UV singularities, 
while the possible IR singularities are still contained in the integral representations 
of the finite parts hi, with the intention of calculating these finite parts by numerical 
integration. 

Therefore it is useful to separate first the infrared part of the two-loop diagrams 
in an analytically manageable form (one-loop diagrams in this case). The "real two- 
loop" calculation which is left after this separation is then free of infrared singularities. 
The analytical separation of infrared divergencies is performed by noticing that the IR 
behaviour of the two-loop diagrams comes essentially only from the loop integration 
over the gluon propagator. Then, the IR behaviour is left unchanged if the loop 
momentum on the propagator common to the two loop integrations (r) is being 
"frozen" to the loop momentum of the IR-finite loop. This analytical separation of 
the IR behaviour is given in figure |] for all IR divergent two-loop diagrams involved 
in this process. 

Once the IR separation is performed, the IR finite part (V2 — V^ 1 etc.) can 
be calculated by numerical integration. For the numerical integration to be stable, 

( J ¥f\ 

one must make sure that the two components, e.g. Vi and V 2 , have the same 
Feynman parameterization such that the singularities cancel already before the Feyn- 
man parameters integration. Then, the IR divergent part (e.g. ) can easily be 

treated analytically, in general by using dimensional regularization to regularize the 
IR singularities. 

For the case of the process considered in this example, a gluon mass regulator can 
also be used. This is possible because in this particular order in a s the IR structure 
is the same as in the Abelian case, and a gluon mass regulator can be used without 
upsetting the Slavnov- Taylor identities. 

It is still an open question if such an analytical separation of the infrared singu- 
larities can always be performed such that the numerical integration can be carried 
out, in the case of other processes with potentially more complicated IR structure. 
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Figure 8: Extracting the infrared divergent pieces of the two-loop diagrams analytically. 
The infrared divergency of the two-loop diagram is the same as the infrared divergency 
of the product of the two one-loop diagrams obtained by 'freezing" the common line 
in the loop momenta integration. 
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This question clearly requires further investigation. 
4.3 Numerical results 

We subtract the UV divergence of the diagrams of figure [7| by performing minimal 
subtractions of overall-divergences and sub-divergences. Further, we extract the IR 
divergences analytically, according to the discussion of the previous section. The ten- 
sor decomposition results in a set of convergent integrals over the remaining Feynman 
parameters, which can be carried out numerically. 

Because in the top mass range of physical interest the diagrams are not close to 
a threshold, the top mass dependency in the relevant range is relatively smooth, as 
it can be seen in figure |9|. We give numerical values for these diagrams in table 1 
for three values of the top mass — intermediary values can be well approximated by 
interpolation. In performing the numerical integration, we took Myy = 80.41 GeV, 
Mz = 91.187 GeV, and the effective electroweak mixing angle siia 2 9w = -232. For 
simplifying both the analytical and the numerical work, we neglected the b mass, and 
in this limit all diagrams become proportional to 7^(1 — 75). It is perfectly possible 
to keep the exact m b dependency throughout the whole calculation at the price of 
dealing with longer expressions and more computing time; however the b mass effect 
on the two-loop Z — > bb decay rate is expected to be very small. 

The numerical results are given in table 1. As for the numerical accuracy and 
efficiency which can be attained, the results of table 1 are for a total of 26 individual 
Feynman diagrams (W and <fr exchange counted separately), each of them evaluated 
for 3 different values of the top mass. To carry out the numerical integration for 
this total of 78 two-loop individual Feynman graph evaluations with an accuracy of 
10 -3 , a total of 100 hours computing time on a 600 MHz Pentium machine was used. 
Higher numerical accuracy can be obtained by simply using more computing power. 
The analytical tensor decomposition part, performed by a FORM computer algebra 
program (we also used a Schoonship version with similar results), takes about one 
hour. 

5 Conclusions 

We described an algorithm for the tensor reduction of massive two-loop diagrams. It 
applies in principle to any massive two-loop graph, and it can be automatized in the 
form of a computer algebra program. The tensor decomposition algorithm results in 
a set of ten special functions hi which are defined in terms of one-dimensional integral 
representations. We described the numerical methods which are used for carrying out 
the remaining integrations. 

By applying the analytical reduction and numerical integration to an important 
three-point example, Z — > bb, we have shown that it can be used in realistic calcu- 
lations, where several internal mass and external momenta scales are involved. This 
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Figure 9: Numerical results for the two-loop three-point function diagrams shown in 
figure 0. The UV divergences and sub- divergences are subtracted by minimal subtrac- 
tion. The IR divergences are subtracted according to figure [|, so that the final result 
for diagrams V2, V4, V&, and V7 are the sum of the finite parts given in the plots plus 
the IR divergent one-loop contributions given in figure 0. There is an overall factor 
0/27^(1 — 75)a s (g 3 /i2 cosdyy). The solid line is the real part, and the dashed line is 
the imaginary part. 
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m t = 165 GeV 


mt = 175 GeV 


m t = 185 GeV 


Vi 




-1.009 ■ 10~ 3 


-7.187- 10~ 4 


-4.057- 10~ 4 


v 2 - 


( T £?\ 

- V 2 { R) 


(-2.873 + ^2.122) • 10~ 3 
1.545 • IO- 3 


(-2.490 + 11.147) • 10~ 3 
2.255 • 10~ 3 


(-2.087 + i.09274) • 10~ 3 
3.034- 10~ 3 


v 3 




v 4 - 


- vi IR) 


(1.215 -z2.481) • 10~ 2 


(1.242-i2.570) • 10~ 2 
2.469 • 10~ 2 


(1.266 -i2.660) -10- 2 
2.861 • 10~ 2 


v 5 




2.107- 10" 2 






(3.089 - i4.257) ■ 10~ 2 


(3.500 - i4.824) ■ 10~ 2 


(3.950 - i5.445) • 10~ 2 


v 7 - 


- vi IR) 


(-.7778 + ^1.281) ■ 10~ 2 


(-.8001 +il.349) ■ 10- 2 
-1.474- lO" 3 


(-.8232 + ^1.420) ■ 10~ 2 






-1.059- 10~ 3 


-1.942 • 10~ 3 


v 9 




6.289 • 10~ 2 


6.703 • 10~ 2 


7.143 • 10- 2 


v l0 




-1.402 • 10~ 2 


-1.389- 10" 2 


-1.380- 10" 2 



Table 1: Numerical values for the two-loop diagrams shown in figure |7|. Vj-Vio 
are the sums of the corresponding W and exchange graphs. An overall color and 
coupling constant factor of i r y fl (l — , j 5 )a s (g 3 /12 cos 0w) is understood. The UV and IR 
divergences are removed as discussed in the text. The numerical integration accuracy 
is 10 -3 . The evaluation of a total of 78 Feynman graph evaluations with this precision 
requires 100 hours computing time on a 600 MHz Pentium machine. 

approach works for any such combination of kinematic variables, apart maybe from 
possible infrared complications. 

In the context of the Z — > bb example given in this paper, we discussed the 
analytical separation of the infrared divergencies. Within our two-loop methods, 
if a process involves infrared singularities, these have to be dealt with in a special 
way because the numerical nature of our methods. It is an open question if this 
IR treatment can be generalized to other two-loop situations with potentially more 
complicated IR structure. 
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Appendix A 

In section 2.2 we introduced two types of scalar functions which are involved in the 
tensor decomposition: Hi and V^ b ia2a3 . Then, we gave integral representations only 
for Hi- The necessary V^w scalar functions can be derived from the corresponding 
"Hi by partial fractioning: 
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where Ti and T 2 are the Euclidian one-loop tadpole integrals: 
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In section 2.4, in order to simplify the discussion, we introduced the functions 
Va b ia2a3 which are related to V^ a2as by a simple loop momentum shift: 
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The functions with ol\ = a 2 = a 3 = 1 which appear in the relations 18 can be cal- 
culated by partial p (eq. 15). For simplifying the notation, we omitted in the above 
formulae the mass and momentum arguments of the functions Ti^rrii, m 2 , m 3 ; k 2 ), 
V^ a2a3 {m 1 ,m 2 ,m 3 , k 2 ), and a2 Q3 ("7.1 , m 2 , m 3 ; k 2 ), and understand that these ar- 
guments appear in this same order in all relations. 

Appendix B 

The integral representations of the ten special functions functions hi, given in eqs. 8, 
are built from the following functions: 



g(m l ,m 2 ,m 3 ;k ;x) = Sp(- ) + Sp(- ) + y x lo. 
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where we use the following notations: 
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and 
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In the above expressions, one special case must be treated separately, namely 
k 2 = 0. The hi functions have a smooth limit for k 2 — > 0. For the purpose of 
numerical evaluation, it is useful to use a Taylor expansion of the functions g, fx, f2, 
and fs around k 2 = for extremely small values of k 2 , where a direct evaluation by 
means of the exact expressions given above would be affected by large cancellations. 
It can be checked that this limit is regular and our approach reduces to the functions 
introduced by van der Bij and Veltman in ref. [jlOfl . 
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